%% 最优化计算方法 - 第一章演示
% 梯度下降法求解无约束优化问题
% 以二次函数f(x) = x1^2 + 2*x2^2为例

% 清空工作区和命令窗口
clear all;
clc;

% 定义目标函数
f = @(x) x(1)^2 + 2*x(2)^2;

% 定义梯度函数
grad_f = @(x) [2*x(1); 4*x(2)];

% 梯度下降法参数
max_iter = 100;      % 最大迭代次数
tol = 1e-6;          % 收敛容差
alpha = 0.1;         % 步长（学习率）

% 初始点
x = [1; 2];
iter = 0;
grad_norm = inf;

% 记录迭代历史
x_history = x;
f_history = f(x);

% 梯度下降迭代
fprintf('迭代开始:\n');
fprintf('迭代次数\t\tf(x)\t\t||∇f(x)||\n');
fprintf('%d\t\t%.6f\t%.6f\n', iter, f(x), norm(grad_f(x)));

while iter < max_iter && grad_norm > tol
    % 计算梯度
    grad = grad_f(x);
    grad_norm = norm(grad);
    
    % 更新位置
    x = x - alpha * grad;
    
    % 更新迭代次数
    iter = iter + 1;
    
    % 记录历史
    x_history = [x_history, x];
    f_history = [f_history, f(x)];
    
    % 打印当前迭代信息
    fprintf('%d\t\t%.6f\t%.6f\n', iter, f(x), grad_norm);
end

% 结果输出
if grad_norm <= tol
    fprintf('\n梯度下降法收敛!\n');
else
    fprintf('\n达到最大迭代次数，未收敛!\n');
end

fprintf('最优解: x1 = %.6f, x2 = %.6f\n', x(1), x(2));
fprintf('最优值: f(x) = %.6f\n', f(x));

% 可视化迭代过程
figure;
[X1, X2] = meshgrid(linspace(-1, 1, 50), linspace(-2, 2, 50));
Z = zeros(size(X1));
for i = 1:size(X1, 1)
    for j = 1:size(X1, 2)
        Z(i, j) = f([X1(i, j); X2(i, j)]);
    end
end

% 绘制等高线
contour(X1, X2, Z, 20);
hold on;

% 绘制迭代路径
plot(x_history(1, :), x_history(2, :), 'r.-', 'MarkerSize', 10);
plot(x_history(1, 1), x_history(2, 1), 'go', 'MarkerSize', 10); % 起点
plot(x_history(1, end), x_history(2, end), 'rs', 'MarkerSize', 10); % 终点

title('梯度下降法优化过程');
xlabel('x_1');
ylabel('x_2');
legend('等高线', '迭代路径', '初始点', '最优点', 'Location', 'best');
grid on;

% 绘制目标函数值随迭代次数的变化
figure;
plot(0:iter, f_history, 'b.-', 'LineWidth', 1.5);
title('目标函数值随迭代次数的变化');
xlabel('迭代次数');
ylabel('f(x)');
grid on; 